function [ newSkelImage ] = fixLoops( skelImage, minLoopSize )
%FIXLOOP Replaces loops in a skeleton with single lines

%Replaces loops with area of the hole below a threshold with lines from all
%branchpoints on the loop to the loop centroid.

    newSkelImage = skelImage;
   
    %Get boundaries of all objects and holes 
    [ boundary,~, n,~] = bwboundaries(skelImage,'holes');
        
    %If there are holes
    if size(boundary,1) > n
        
        %For each hole (holes are in positions greater than n in the cell
        %araay
        for i = n+1:length(boundary)
             theHoleBoundary = boundary{i};
             
             %Create a mask of the hole boundary
             mask = zeros(size(skelImage));
             mask( sub2ind(size(mask), theHoleBoundary(:,1), theHoleBoundary(:,2) )) = 1;
              
             %Get the area of the hole
             mask = imfill(mask);
             areaOfHole = length(find(mask));
                 
             %If the hole is smaller than the min loop size, then replace
             %it with lines from the branchpoints to the centroid
                if areaOfHole < minLoopSize
                    
                    %Thicken the mask of the hole by one pixel so it covers the
                    %loop
                    nhood =[1 1 1; 1 1 1; 1 1 1];
                    mask = imdilate(mask,nhood, 'same');
                               
                    %The mask applied to the skel image will give the loop
                    %pixel image.
                    loopImage = mask & skelImage;
                    loopRegion = regionprops(loopImage, 'PixelList','PixelIdxList', 'Area','Centroid');

                
                    %Find the branchpoints in the original skel Image that
                    %are also part of the loop
                    branchImage = bwmorph(skelImage, 'branchpoints');
                    imBranchesOnLoop = branchImage&loopImage;
                    imNonBranchLoopPixels = loopImage & (~ imBranchesOnLoop);
                
                    %Get the row and column of branch pixels on the loop
                    [row, col] = find(imBranchesOnLoop);
                    branchPixel = [row, col];
                    
                    %Determine center of loop
                    centroidRow = round(loopRegion.Centroid(2));
                    centroidCol = round(loopRegion.Centroid(1));
                    centerOfLoop = [centroidRow, centroidCol];
                 
                               
                    %Get indices of all lines from all branch points on the
                    %loop to the centroid.
                    newLineIndices = [];
                    for j = 1:size(branchPixel,1);

                        index = getLineSegmentIndices(skelImage,...
                                       centerOfLoop, branchPixel);
                        newLineIndices = [newLineIndices, index];
                    end
                
                    %Delete non branch pixels from the original image and add
                    %in lines from each branch to the loop centroid.
                    newSkelImage = newSkelImage - imNonBranchLoopPixels;
                    newSkelImage(newLineIndices) = 1;
                    newSkelImage(centroidRow, centroidCol) = 1;
                
                end
             
           
        end
    end
             
    

end

